
*********************************************************************************
** Run models and create graphs in all non-metro counties 
*********************************************************************************
set more off
use "Huntingseason_Analysisfile_Dataverse", clear

** Models for all shootings
poisson nshot ib3.hwkstart ib2014.year ib1.dow ib10.month tg xm, cluster(countyid) irr
estimates store nshotnofe2
xtpoisson nshot ib3.hwkstart ib2014.year ib1.dow ib10.month tg xm, fe i(countyid) vce(robust) irr
estimates store nshot2
xtpoisson nshotnh ib3.hwkstart ib2014.year ib1.dow ib10.month tg xm, fe i(countyid) vce(robust) irr
estimates store nshotnh2

set scheme s2color

** Plot without notes
coefplot ///
	(nshotnofe2, eform keep(*.hwkstart) base ciopts(lcolor(gs0)) ///
	yline(1, lpattern(dash) lcolor(gs10)) xline(3, lpattern(dash) lcolor(gs10)) ///
	graphregion(color(white)) bgcolor(white)  ///
	xscale(lcolor(gs3)) yscale(log lcolor(gs3) range(.5 2)) ///
	xlabel(, labcolor("gs1")) ylabel(.5 .67 1 1.5 2, labcolor("gs1"))  ///
	msymbol(d) msize(small) mlcolor(gs0) mfcolor(gs0) label("No fixed effects")) ///
	(nshot2, eform ///
	keep(*.hwkstart) base ciopts(lcolor(gs6)) ///
	msymbol(s) msize(small) mlcolor(gs6) mfcolor(gs6) label("County fixed effects")) ///
	(nshotnh2, eform ///
	keep(*.hwkstart) base ciopts(lcolor(gs10)) ///
	msymbol(circle) msize(small) mlcolor(gs10) mfcolor(gs10) label("Exclude hunting accidents")) ///
	, legend(size(2) rows(1) position(12)) vertical ytitle("Incidence rate ratio, log scale, log scale", size(small)) ylabel(, labsize(vsmall)) ///
	xlabel(1 "3 weeks before" 2 "2 weeks before" 3 "1 week before" 4 "week 1 deer season" 5 "week 2" 6 "week 3", labsize(vsmall) angle(vertical)) generate(Fig2_)
	
** Check whether relationship is stronger for hand guns or long guns 
xtpoisson nshotlg ib3.hwkstart ib2014.year ib1.dow ib10.month tg xm, fe i(countyid) vce(robust) irr
xtpoisson nshothg ib3.hwkstart ib2014.year ib1.dow ib10.month tg xm, fe i(countyid) vce(robust) irr

** Run models for fatal shootings
poisson nkilled ib3.hwkstart ib2014.year ib1.dow ib10.month tg xm, cluster(countyid) irr
estimates store nkillednofe2
xtpoisson nkilled ib3.hwkstart ib2014.year ib1.dow ib10.month tg xm, fe i(countyid) vce(robust) irr
estimates store nkilled2
xtpoisson nkillednh ib3.hwkstart ib2014.year ib1.dow ib10.month tg xm, fe i(countyid) vce(robust) irr
estimates store nkillednh2

** Plot 
coefplot ///
	(nkillednofe2, eform keep(*.hwkstart) base ciopts(lcolor(gs0)) ///
	yline(1, lpattern(dash) lcolor(gs10)) xline(3, lpattern(dash) lcolor(gs10)) ///
	graphregion(color(white)) bgcolor(white)  ///
	xscale(lcolor(gs3)) yscale(log lcolor(gs3) range(.5 2)) ///
	xlabel(, labcolor("gs1")) ylabel(.5 .67 1 1.5 2, labcolor("gs1"))  ///
	msymbol(d) msize(small) mlcolor(gs0) mfcolor(gs0) label("No fixed effects")) ///
	(nkilled2, eform ///
	keep(*.hwkstart) base ciopts(lcolor(gs6)) ///
	msymbol(s) msize(small) mlcolor(gs6) mfcolor(gs6) label("County fixed effects")) ///
	(nkillednh2, eform ///
	keep(*.hwkstart) base ciopts(lcolor(gs10)) ///
	msymbol(circle) msize(small) mlcolor(gs10) mfcolor(gs10) label("Exclude hunting accidents")) ///
	, legend(size(2) rows(1) position(12)) vertical ytitle("Incidence rate ratio, log scale", size(small)) ylabel(, labsize(vsmall)) ///
	xlabel(1 "3 weeks before" 2 "2 weeks before" 3 "1 week before" 4 "week 1 deer season" 5 "week 2" 6 "week 3", labsize(vsmall) angle(vertical))
	
xtpoisson nshot ib3.hwkstart ib2014.year ib1.dow ib10.month tg xm if mhuntlic3==1, fe i(countyid) vce(robust) irr
estimates store nshot1b
xtpoisson nshot ib3.hwkstart ib2014.year ib1.dow ib10.month tg xm if mhuntlic3==2, fe i(countyid) vce(robust) irr
estimates store nshot2b
xtpoisson nshot ib3.hwkstart ib2014.year ib1.dow ib10.month tg xm if mhuntlic3==3, fe i(countyid) vce(robust) irr
estimates store nshot3b

** Plot
coefplot ///
	(nshot1b, eform keep(*.hwkstart) base ciopts(lcolor(gs0)) ///
	yline(1, lpattern(dash) lcolor(gs10)) xline(3, lpattern(dash) lcolor(gs10)) ///
	graphregion(color(white)) bgcolor(white)  ///
	xscale(lcolor(gs3)) yscale(log lcolor(gs3)) ///
	xlabel(, labcolor("gs1")) ylabel(.5 1 2 4, labcolor("gs1"))  ///
	msymbol(d) msize(small) mlcolor(gs0) mfcolor(gs0) label("Bottom third hunting licenses per capita")) ///
	(nshot2b, eform ///
	keep(*.hwkstart) base ciopts(lcolor(gs6)) ///
	msymbol(s) msize(small) mlcolor(gs6) mfcolor(gs6) label("Middle third hunting licenses per capita")) ///
	(nshot3b, eform ///
	keep(*.hwkstart) base ciopts(lcolor(gs10)) ///
	msymbol(circle) msize(small) mlcolor(gs10) mfcolor(gs10) label("Top third hunting licenses per capita")) ///
	, legend(size(2) rows(1) position(12)) vertical ytitle("Incidence rate ratio, log scale", size(small)) ylabel(, labsize(vsmall)) ///
	xlabel(1 "3 weeks before" 2 "2 weeks before" 3 "1 week before" 4 "week 1 deer season" 5 "week 2" 6 "week 3", labsize(vsmall) angle(vertical)) generate(Fig3_)
